Jan06AssignmentQs
Xiaoyan Yue
1/5/2017
Exercises
To gain familiarity with brms I am going to have you refit models from Statistical Rethinking with brms. In all cases compare the brms and map2stan parameter estimates; they should be effectively the same if you have specified them correctly.
Q1:Fit model 10.9 (R code 10.28) with brms. Remember that you should not need to transform the predictors, create dummy variables, nor coerce indexes to use brms. Compare to the map2stan fit. Test whether the coefficient for gender is different from 0 in the brms model.
#Load the data
data(UCBadmit)
d <- UCBadmit
d## dept applicant.gender admit reject applications
## 1 A male 512 313 825
## 2 A female 89 19 108
## 3 B male 353 207 560
## 4 B female 17 8 25
## 5 C male 120 205 325
## 6 C female 202 391 593
## 7 D male 138 279 417
## 8 D female 131 244 375
## 9 E male 53 138 191
## 10 E female 94 299 393
## 11 F male 22 351 373
## 12 F female 24 317 341summary(d)## dept applicant.gender admit reject applications
## A:2 female:6 Min. : 17.00 Min. : 8.0 Min. : 25.0
## B:2 male :6 1st Qu.: 45.75 1st Qu.:188.2 1st Qu.:291.5
## C:2 Median :107.00 Median :261.5 Median :374.0
## D:2 Mean :146.25 Mean :230.9 Mean :377.2
## E:2 3rd Qu.:154.00 3rd Qu.:314.0 3rd Qu.:452.8
## F:2 Max. :512.00 Max. :391.0 Max. :825.0#Fit the model with map2stan
#First, creat a dummy variable for "application.gender"
d$male <- ifelse( d$applicant.gender=="male" , 1 , 0 )
#Second, make index for each dpartment
d$dept_id <- coerce_index( d$dept )
#Third, get rid of "." in the colum name "application.gender"
colnames(d) <- sub(".","_",colnames(d),fixed = TRUE)
mQ1.stan <- map2stan(
alist(
admit ~ dbinom( applications , p ) ,
logit(p) <- a_dept[dept_id] + bm*male ,
a_dept[dept_id] ~ dnorm(0,10) ,
bm ~ dnorm(0,10)
) , data=d, chains = 4)## In file included from file10a4d398053.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math/rev/core.hpp:12:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math/rev/core/gevv_vvv_vari.hpp:5:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math/rev/core/var.hpp:7:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/BH/include/boost/math/tools/config.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/BH/include/boost/config.hpp:39:
## /Library/Frameworks/R.framework/Versions/3.3/Resources/library/BH/include/boost/config/compiler/clang.hpp:196:11: warning: 'BOOST_NO_CXX11_RVALUE_REFERENCES' macro redefined [-Wmacro-redefined]
## # define BOOST_NO_CXX11_RVALUE_REFERENCES
## ^
## <command line>:6:9: note: previous definition is here
## #define BOOST_NO_CXX11_RVALUE_REFERENCES 1
## ^
## 1 warning generated.
##
## SAMPLING FOR MODEL 'admit ~ dbinom(applications, p)' NOW (CHAIN 1).
##
## Chain 1, Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1, Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1, Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1, Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1, Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1, Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1, Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1, Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1, Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1, Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1, Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1, Iteration: 2000 / 2000 [100%] (Sampling)
## Elapsed Time: 0.050203 seconds (Warm-up)
## 0.038984 seconds (Sampling)
## 0.089187 seconds (Total)
##
##
## SAMPLING FOR MODEL 'admit ~ dbinom(applications, p)' NOW (CHAIN 2).
##
## Chain 2, Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2, Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2, Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2, Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2, Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2, Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2, Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2, Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2, Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2, Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2, Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2, Iteration: 2000 / 2000 [100%] (Sampling)
## Elapsed Time: 0.048116 seconds (Warm-up)
## 0.044374 seconds (Sampling)
## 0.09249 seconds (Total)
##
##
## SAMPLING FOR MODEL 'admit ~ dbinom(applications, p)' NOW (CHAIN 3).
##
## Chain 3, Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3, Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3, Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3, Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3, Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3, Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3, Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3, Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3, Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3, Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3, Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3, Iteration: 2000 / 2000 [100%] (Sampling)
## Elapsed Time: 0.051587 seconds (Warm-up)
## 0.05759 seconds (Sampling)
## 0.109177 seconds (Total)
##
##
## SAMPLING FOR MODEL 'admit ~ dbinom(applications, p)' NOW (CHAIN 4).
##
## Chain 4, Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4, Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4, Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4, Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4, Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4, Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4, Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4, Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4, Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4, Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4, Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4, Iteration: 2000 / 2000 [100%] (Sampling)
## Elapsed Time: 0.056264 seconds (Warm-up)
## 0.04476 seconds (Sampling)
## 0.101024 seconds (Total)
##
##
## SAMPLING FOR MODEL 'admit ~ dbinom(applications, p)' NOW (CHAIN 1).
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
##
## Chain 1, Iteration: 1 / 1 [100%] (Sampling)
## Elapsed Time: 3e-06 seconds (Warm-up)
## 3.8e-05 seconds (Sampling)
## 4.1e-05 seconds (Total)## Computing WAIC## Constructing posterior predictions## [ 400 / 4000 ]
[ 800 / 4000 ]
[ 1200 / 4000 ]
[ 1600 / 4000 ]
[ 2000 / 4000 ]
[ 2400 / 4000 ]
[ 2800 / 4000 ]
[ 3200 / 4000 ]
[ 3600 / 4000 ]
[ 4000 / 4000 ]## Aggregated binomial counts detected. Splitting to 0/1 outcome for WAIC calculation.#test the model
precis(mQ1.stan,depth = 2)## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## a_dept[1] 0.68 0.10 0.51 0.83 2799 1
## a_dept[2] 0.64 0.12 0.46 0.83 3026 1
## a_dept[3] -0.58 0.07 -0.69 -0.46 4000 1
## a_dept[4] -0.61 0.09 -0.74 -0.47 4000 1
## a_dept[5] -1.06 0.10 -1.22 -0.91 4000 1
## a_dept[6] -2.64 0.16 -2.89 -2.38 4000 1
## bm -0.10 0.08 -0.23 0.03 2322 1pairs(mQ1.stan)
summary(mQ1.stan)## Inference for Stan model: admit ~ dbinom(applications, p).
## 4 chains, each with iter=2000; warmup=1000; thin=1;
## post-warmup draws per chain=1000, total post-warmup draws=4000.
##
## mean se_mean sd 2.5% 25% 50% 75%
## a_dept[1] 0.68 0.00 0.10 0.49 0.62 0.68 0.75
## a_dept[2] 0.64 0.00 0.12 0.41 0.56 0.64 0.72
## a_dept[3] -0.58 0.00 0.07 -0.73 -0.63 -0.58 -0.53
## a_dept[4] -0.61 0.00 0.09 -0.79 -0.67 -0.61 -0.56
## a_dept[5] -1.06 0.00 0.10 -1.26 -1.13 -1.06 -0.99
## a_dept[6] -2.64 0.00 0.16 -2.96 -2.74 -2.63 -2.52
## bm -0.10 0.00 0.08 -0.26 -0.16 -0.10 -0.05
## dev 96.13 0.09 3.80 90.78 93.25 95.48 98.26
## lp__ -2597.29 0.04 1.90 -2601.92 -2598.35 -2596.96 -2595.85
## 97.5% n_eff Rhat
## a_dept[1] 0.89 2799 1
## a_dept[2] 0.87 3026 1
## a_dept[3] -0.44 4000 1
## a_dept[4] -0.44 4000 1
## a_dept[5] -0.87 4000 1
## a_dept[6] -2.34 4000 1
## bm 0.05 2322 1
## dev 105.38 1917 1
## lp__ -2594.61 1917 1
##
## Samples were drawn using NUTS(diag_e) at Thu Jan 5 18:19:29 2017.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).#Fit the model by brms, I'm confused by how to create the outcome variable, so I tried different ways
#(in chapter10, 10.1.4, an aggregated binomial uses cbind to build the outcome variable)
mQ1.brms1 <- brm(cbind(admit,reject) ~ 0 + dept + applicant_gender,#separate intercept for each dept
#family = "binomial",#Multivariate models are not yet implemented for family 'binomial'.
prior=set_prior("normal(0,10)", class="b"),
data = d)## Compiling the C++ model## Start samplingsummary(mQ1.brms1)## Family: gaussian (identity)
## Formula: cbind(admit, reject) ~ 0 + dept + applicant_gender
## Data: d (Number of observations: 12)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
## WAIC: Not computed
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI
## admit_deptA 2.11 10.04 -17.51 21.77
## admit_deptB 1.32 10.06 -18.32 21.73
## admit_deptC 0.04 9.86 -18.69 19.58
## admit_deptD 0.08 10.10 -19.49 20.17
## admit_deptE -0.38 10.09 -19.65 19.14
## admit_deptF -1.38 9.62 -20.19 17.57
## admit_applicant_gendermale 2.14 9.96 -17.45 21.62
## reject_deptA -0.43 9.73 -19.59 18.84
## reject_deptB -0.26 9.65 -19.48 18.51
## reject_deptC 0.99 10.05 -18.69 20.78
## reject_deptD 0.89 9.77 -19.19 20.29
## reject_deptE 1.03 10.02 -18.70 21.04
## reject_deptF 1.89 9.83 -16.94 21.22
## reject_applicant_gendermale 1.47 10.12 -18.04 21.42
## Eff.Sample Rhat
## admit_deptA 4000 1
## admit_deptB 4000 1
## admit_deptC 4000 1
## admit_deptD 4000 1
## admit_deptE 4000 1
## admit_deptF 4000 1
## admit_applicant_gendermale 4000 1
## reject_deptA 4000 1
## reject_deptB 4000 1
## reject_deptC 4000 1
## reject_deptD 4000 1
## reject_deptE 4000 1
## reject_deptF 4000 1
## reject_applicant_gendermale 4000 1
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sigma(admit) 203.77 39.88 143.53 298.03 4000 1
## sigma(reject) 256.81 52.83 178.56 377.04 4000 1
## rescor(admit,reject) 0.61 0.17 0.19 0.86 4000 1
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).plot(mQ1.brms1, ask = FALSE)
#then, I create a outcome variable for the model, calculate the probablity of admit in each department
d$admit_p <- d$admit/d$applications
mQ1.brms2 <- brm(admit_p ~ 0 + dept + applicant_gender,#separate intercept for each dept
prior=set_prior("normal(0,10)", class="b"),
data = d)## Compiling the C++ model
## Start samplingsummary(mQ1.brms2)## Family: gaussian (identity)
## Formula: admit_p ~ 0 + dept + applicant_gender
## Data: d (Number of observations: 12)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
## WAIC: Not computed
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## deptA 0.74 0.09 0.60 0.89 4000 1
## deptB 0.67 0.08 0.53 0.82 3748 1
## deptC 0.37 0.08 0.23 0.52 3762 1
## deptD 0.36 0.08 0.20 0.50 3424 1
## deptE 0.28 0.08 0.13 0.43 3135 1
## deptF 0.08 0.08 -0.07 0.24 2632 1
## applicant_gendermale -0.04 0.06 -0.14 0.08 2942 1
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sigma 0.09 0.06 0.04 0.21 483 1.01
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).plot(mQ1.brms2, ask = FALSE)
#Test whether the coefficient for gender is different from 0 in the brms model
hypothesis(mQ1.brms1,"admit_applicant_gendermale = 0") ## Hypothesis Tests for class b:
## Estimate Est.Error l-95% CI u-95% CI Evid.Ratio
## (admit_applicant_... = 0 2.14 9.96 -17.45 21.62 NA
## ---
## '*': The expected value under the hypothesis lies outside the 95% CI.plot(hypothesis(mQ1.brms1,"admit_applicant_gendermale = 0"))
hypothesis(mQ1.brms2,"applicant_gendermale = 0")## Hypothesis Tests for class b:
## Estimate Est.Error l-95% CI u-95% CI Evid.Ratio
## (applicant_gender... = 0 -0.04 0.06 -0.14 0.08 NA
## ---
## '*': The expected value under the hypothesis lies outside the 95% CI.plot(hypothesis(mQ1.brms2,"applicant_gendermale = 0"))
#fit model with interactions between department and gender
#(1)
mQ1.brms3 <- brm(cbind(admit,reject) ~ 0 + dept * applicant_gender,
prior=set_prior("normal(0,10)", class="b"),
data = d)## Compiling the C++ model
## Start samplingsummary(mQ1.brms3)## Family: gaussian (identity)
## Formula: cbind(admit, reject) ~ 0 + dept * applicant_gender
## Data: d (Number of observations: 12)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
## WAIC: Not computed
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI
## admit_deptA 2.06 10.15 -17.47 21.87
## admit_deptB 1.23 10.16 -18.30 20.62
## admit_deptC -0.10 10.09 -19.59 19.70
## admit_deptD -0.09 10.03 -19.64 19.18
## admit_deptE -0.25 10.18 -19.82 19.41
## admit_deptF -1.31 9.96 -20.82 18.62
## admit_applicant_gendermale 2.19 9.69 -16.89 20.69
## admit_deptB:applicant_gendermale 1.06 9.94 -18.19 20.13
## admit_deptC:applicant_gendermale -0.06 9.96 -18.90 19.25
## admit_deptD:applicant_gendermale 0.08 9.56 -18.49 18.82
## admit_deptE:applicant_gendermale 0.02 10.16 -19.75 19.80
## admit_deptF:applicant_gendermale -0.70 9.89 -19.97 18.17
## reject_deptA -0.61 9.62 -19.22 17.86
## reject_deptB -0.36 9.88 -19.93 18.91
## reject_deptC 1.07 10.08 -18.49 20.72
## reject_deptD 1.01 9.82 -18.33 20.24
## reject_deptE 0.86 9.93 -18.64 20.32
## reject_deptF 1.79 9.79 -17.66 20.80
## reject_applicant_gendermale 1.62 9.88 -17.36 20.90
## reject_deptB:applicant_gendermale -0.05 9.83 -19.20 19.27
## reject_deptC:applicant_gendermale 0.23 10.08 -19.14 20.00
## reject_deptD:applicant_gendermale 0.54 9.54 -18.35 19.03
## reject_deptE:applicant_gendermale 0.42 9.51 -18.18 19.38
## reject_deptF:applicant_gendermale 1.01 10.00 -18.94 20.76
## Eff.Sample Rhat
## admit_deptA 4000 1
## admit_deptB 4000 1
## admit_deptC 4000 1
## admit_deptD 4000 1
## admit_deptE 4000 1
## admit_deptF 4000 1
## admit_applicant_gendermale 4000 1
## admit_deptB:applicant_gendermale 4000 1
## admit_deptC:applicant_gendermale 4000 1
## admit_deptD:applicant_gendermale 4000 1
## admit_deptE:applicant_gendermale 4000 1
## admit_deptF:applicant_gendermale 4000 1
## reject_deptA 4000 1
## reject_deptB 4000 1
## reject_deptC 4000 1
## reject_deptD 4000 1
## reject_deptE 4000 1
## reject_deptF 4000 1
## reject_applicant_gendermale 4000 1
## reject_deptB:applicant_gendermale 4000 1
## reject_deptC:applicant_gendermale 4000 1
## reject_deptD:applicant_gendermale 4000 1
## reject_deptE:applicant_gendermale 4000 1
## reject_deptF:applicant_gendermale 4000 1
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sigma(admit) 203.69 41.01 142.10 303.57 4000 1
## sigma(reject) 256.18 50.63 179.50 372.64 4000 1
## rescor(admit,reject) 0.61 0.17 0.21 0.85 4000 1
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).plot(mQ1.brms3, ask = FALSE)
#(2)
mQ1.brms4 <- brm(admit_p ~ 0 + dept * applicant_gender,#separate intercept for each dept
prior=set_prior("normal(0,10)", class="b"),
data = d)## Compiling the C++ model
## Start sampling## Warning: There were 402 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup## Warning: There were 4 chains where the estimated Bayesian Fraction of Missing Information was low. See
## http://mc-stan.org/misc/warnings.html#bfmi-low## Warning: Examine the pairs() plot to diagnose sampling problemssummary(mQ1.brms4)## Warning: The model has not converged (some Rhats are > 1.1). Do not analyse the results!
## We recommend running more iterations and/or setting stronger priors.## Family: gaussian (identity)
## Formula: admit_p ~ 0 + dept * applicant_gender
## Data: d (Number of observations: 12)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
## WAIC: Not computed
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample
## deptA 0.72 1.93 -3.75 4.79 2615
## deptB 0.69 2.14 -3.78 5.64 2925
## deptC 0.28 2.14 -4.33 4.74 3531
## deptD 0.40 2.25 -4.45 5.44 3637
## deptE 0.24 2.13 -4.45 4.93 3235
## deptF 0.13 2.16 -4.59 4.54 4000
## applicant_gendermale -0.08 2.19 -5.00 4.74 1754
## deptB:applicant_gendermale -0.03 3.38 -7.29 7.24 2246
## deptC:applicant_gendermale 0.16 3.37 -7.10 7.34 2234
## deptD:applicant_gendermale -0.04 3.46 -7.48 7.75 2322
## deptE:applicant_gendermale 0.12 3.29 -6.89 7.55 2200
## deptF:applicant_gendermale 0.05 3.60 -7.62 8.25 2106
## Rhat
## deptA 1
## deptB 1
## deptC 1
## deptD 1
## deptE 1
## deptF 1
## applicant_gendermale 1
## deptB:applicant_gendermale 1
## deptC:applicant_gendermale 1
## deptD:applicant_gendermale 1
## deptE:applicant_gendermale 1
## deptF:applicant_gendermale 1
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sigma 1.89 1.57 0.21 5.71 22 1.16
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).plot(mQ1.brms4, ask = FALSE)
#Test whether the coefficient for gender is different from 0 in the brms model
hypothesis(mQ1.brms3,"admit_applicant_gendermale = 0") ## Hypothesis Tests for class b:
## Estimate Est.Error l-95% CI u-95% CI Evid.Ratio
## (admit_applicant_... = 0 2.19 9.69 -16.89 20.69 NA
## ---
## '*': The expected value under the hypothesis lies outside the 95% CI.plot(hypothesis(mQ1.brms3,"admit_applicant_gendermale = 0"))
hypothesis(mQ1.brms4,"applicant_gendermale = 0")## Hypothesis Tests for class b:
## Estimate Est.Error l-95% CI u-95% CI Evid.Ratio
## (applicant_gender... = 0 -0.08 2.19 -5 4.74 NA
## ---
## '*': The expected value under the hypothesis lies outside the 95% CI.plot(hypothesis(mQ1.brms4,"applicant_gendermale = 0"))
###Q2:Refit models 12.1 and 12.2 (Rcode 12.2 and 12.3) with brms
data(reedfrogs)
d2 <- reedfrogs
d2## density pred size surv propsurv
## 1 10 no big 9 0.9000000
## 2 10 no big 10 1.0000000
## 3 10 no big 7 0.7000000
## 4 10 no big 10 1.0000000
## 5 10 no small 9 0.9000000
## 6 10 no small 9 0.9000000
## 7 10 no small 10 1.0000000
## 8 10 no small 9 0.9000000
## 9 10 pred big 4 0.4000000
## 10 10 pred big 9 0.9000000
## 11 10 pred big 7 0.7000000
## 12 10 pred big 6 0.6000000
## 13 10 pred small 7 0.7000000
## 14 10 pred small 5 0.5000000
## 15 10 pred small 9 0.9000000
## 16 10 pred small 9 0.9000000
## 17 25 no big 24 0.9600000
## 18 25 no big 23 0.9200000
## 19 25 no big 22 0.8800000
## 20 25 no big 25 1.0000000
## 21 25 no small 23 0.9200000
## 22 25 no small 23 0.9200000
## 23 25 no small 23 0.9200000
## 24 25 no small 21 0.8400000
## 25 25 pred big 6 0.2400000
## 26 25 pred big 13 0.5200000
## 27 25 pred big 4 0.1600000
## 28 25 pred big 9 0.3600000
## 29 25 pred small 13 0.5200000
## 30 25 pred small 20 0.8000000
## 31 25 pred small 8 0.3200000
## 32 25 pred small 10 0.4000000
## 33 35 no big 34 0.9714286
## 34 35 no big 33 0.9428571
## 35 35 no big 33 0.9428571
## 36 35 no big 31 0.8857143
## 37 35 no small 31 0.8857143
## 38 35 no small 35 1.0000000
## 39 35 no small 33 0.9428571
## 40 35 no small 32 0.9142857
## 41 35 pred big 4 0.1142857
## 42 35 pred big 12 0.3428571
## 43 35 pred big 13 0.3714286
## 44 35 pred big 14 0.4000000
## 45 35 pred small 22 0.6285714
## 46 35 pred small 12 0.3428571
## 47 35 pred small 31 0.8857143
## 48 35 pred small 17 0.4857143# make the tank cluster variable
d2$tank <- 1:nrow(d2)
#fit the model with stan
m12.1.stan <- map2stan(
alist(
surv ~ dbinom( density , p ) ,
logit(p) <- a_tank[tank] ,
a_tank[tank] ~ dnorm( 0 , 5 )
), data=d2 )## In file included from file10a46fa2057e.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math/rev/core.hpp:12:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math/rev/core/gevv_vvv_vari.hpp:5:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math/rev/core/var.hpp:7:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/BH/include/boost/math/tools/config.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/BH/include/boost/config.hpp:39:
## /Library/Frameworks/R.framework/Versions/3.3/Resources/library/BH/include/boost/config/compiler/clang.hpp:196:11: warning: 'BOOST_NO_CXX11_RVALUE_REFERENCES' macro redefined [-Wmacro-redefined]
## # define BOOST_NO_CXX11_RVALUE_REFERENCES
## ^
## <command line>:6:9: note: previous definition is here
## #define BOOST_NO_CXX11_RVALUE_REFERENCES 1
## ^
## 1 warning generated.
##
## SAMPLING FOR MODEL 'surv ~ dbinom(density, p)' NOW (CHAIN 1).
##
## Chain 1, Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1, Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1, Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1, Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1, Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1, Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1, Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1, Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1, Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1, Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1, Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1, Iteration: 2000 / 2000 [100%] (Sampling)
## Elapsed Time: 0.174301 seconds (Warm-up)
## 0.154378 seconds (Sampling)
## 0.328679 seconds (Total)
##
##
## SAMPLING FOR MODEL 'surv ~ dbinom(density, p)' NOW (CHAIN 1).
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
##
## Chain 1, Iteration: 1 / 1 [100%] (Sampling)
## Elapsed Time: 3e-06 seconds (Warm-up)
## 0.0001 seconds (Sampling)
## 0.000103 seconds (Total)## Computing WAIC## Constructing posterior predictions## [ 100 / 1000 ]
[ 200 / 1000 ]
[ 300 / 1000 ]
[ 400 / 1000 ]
[ 500 / 1000 ]
[ 600 / 1000 ]
[ 700 / 1000 ]
[ 800 / 1000 ]
[ 900 / 1000 ]
[ 1000 / 1000 ]## Aggregated binomial counts detected. Splitting to 0/1 outcome for WAIC calculation.precis(m12.1.stan,depth = 2)## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## a_tank[1] 2.54 1.26 0.73 4.36 484 1.00
## a_tank[2] 5.61 2.64 1.62 9.53 586 1.00
## a_tank[3] 0.89 0.72 -0.19 2.00 771 1.00
## a_tank[4] 5.73 2.79 1.63 9.75 666 1.00
## a_tank[5] 2.54 1.23 0.74 4.42 764 1.00
## a_tank[6] 2.58 1.28 0.76 4.50 639 1.00
## a_tank[7] 5.78 2.87 1.85 9.83 407 1.01
## a_tank[8] 2.59 1.16 0.78 4.27 728 1.00
## a_tank[9] -0.45 0.70 -1.61 0.62 1000 1.00
## a_tank[10] 2.57 1.21 0.69 4.37 572 1.00
## a_tank[11] 0.94 0.74 -0.31 2.01 818 1.00
## a_tank[12] 0.44 0.69 -0.60 1.54 1000 1.00
## a_tank[13] 0.93 0.75 -0.34 1.95 1000 1.00
## a_tank[14] 0.01 0.68 -1.08 1.05 1000 1.00
## a_tank[15] 2.48 1.10 0.75 4.06 919 1.00
## a_tank[16] 2.52 1.18 0.80 4.37 481 1.00
## a_tank[17] 3.53 1.21 1.70 5.29 559 1.00
## a_tank[18] 2.56 0.75 1.46 3.73 1000 1.00
## a_tank[19] 2.09 0.65 1.07 3.06 1000 1.00
## a_tank[20] 6.54 2.81 2.50 10.28 510 1.00
## a_tank[21] 2.62 0.81 1.48 3.96 887 1.00
## a_tank[22] 2.64 0.77 1.42 3.78 720 1.00
## a_tank[23] 2.59 0.77 1.51 3.91 785 1.00
## a_tank[24] 1.75 0.59 0.87 2.64 814 1.00
## a_tank[25] -1.19 0.47 -1.91 -0.44 974 1.00
## a_tank[26] 0.10 0.40 -0.52 0.75 1000 1.00
## a_tank[27] -1.70 0.52 -2.51 -0.86 1000 1.00
## a_tank[28] -0.60 0.41 -1.25 0.03 1000 1.00
## a_tank[29] 0.07 0.42 -0.56 0.74 645 1.00
## a_tank[30] 1.41 0.51 0.62 2.18 1000 1.00
## a_tank[31] -0.78 0.45 -1.53 -0.11 920 1.00
## a_tank[32] -0.43 0.41 -1.02 0.32 1000 1.00
## a_tank[33] 3.77 1.08 2.00 5.21 604 1.00
## a_tank[34] 2.94 0.76 1.78 4.08 875 1.00
## a_tank[35] 2.98 0.77 1.80 4.16 866 1.00
## a_tank[36] 2.14 0.56 1.26 2.97 794 1.00
## a_tank[37] 2.13 0.54 1.29 2.96 865 1.00
## a_tank[38] 6.78 2.48 3.13 10.32 652 1.00
## a_tank[39] 2.95 0.80 1.60 4.05 986 1.00
## a_tank[40] 2.48 0.59 1.47 3.31 1000 1.00
## a_tank[41] -2.11 0.52 -2.97 -1.33 808 1.00
## a_tank[42] -0.67 0.36 -1.21 -0.08 1000 1.00
## a_tank[43] -0.54 0.36 -1.07 0.10 1000 1.00
## a_tank[44] -0.43 0.34 -0.97 0.12 885 1.00
## a_tank[45] 0.54 0.35 -0.03 1.08 1000 1.00
## a_tank[46] -0.68 0.36 -1.22 -0.10 672 1.00
## a_tank[47] 2.16 0.57 1.26 3.06 794 1.00
## a_tank[48] -0.06 0.32 -0.59 0.43 1000 1.00#refit the model 12.1 with brms
#(1)without interaction
m12.1.brms1 <- brm(propsurv ~ pred + size,
prior=set_prior("normal(0,5)", class="b"),
data = d2)## Compiling the C++ model## Start samplingsummary(m12.1.brms1)## Family: gaussian (identity)
## Formula: propsurv ~ pred + size
## Data: d2 (Number of observations: 48)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
## WAIC: Not computed
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## Intercept 0.88 0.04 0.79 0.96 4000 1
## predpred -0.40 0.05 -0.50 -0.31 4000 1
## sizesmall 0.09 0.05 -0.01 0.19 4000 1
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sigma 0.17 0.02 0.14 0.22 4000 1
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).plot(m12.1.brms1,ask=FALSE)
#(2)with
m12.1.brms2 <- brm(propsurv ~ pred*size,
prior=set_prior("normal(0,5)", class="b"),
data = d2)## Compiling the C++ model
## Start samplingsummary(m12.1.brms2)## Family: gaussian (identity)
## Formula: propsurv ~ pred * size
## Data: d2 (Number of observations: 48)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
## WAIC: Not computed
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## Intercept 0.93 0.05 0.83 1.02 2793 1
## predpred -0.50 0.07 -0.64 -0.37 2522 1
## sizesmall -0.01 0.07 -0.15 0.13 2449 1
## predpred:sizesmall 0.20 0.10 0.01 0.40 2250 1
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sigma 0.17 0.02 0.14 0.21 3671 1
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).plot(m12.1.brms2,ask=FALSE)
#fit model 12.2 with stan
m12.2.stan <- map2stan(
alist(
surv ~ dbinom( density , p ) ,
logit(p) <- a_tank[tank] ,
a_tank[tank] ~ dnorm( a , sigma ) ,
a ~ dnorm(0,1) ,
sigma ~ dcauchy(0,1)
), data=d2 , chains=4 )## In file included from file10a478873559.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math/rev/core.hpp:12:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math/rev/core/gevv_vvv_vari.hpp:5:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math/rev/core/var.hpp:7:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/BH/include/boost/math/tools/config.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/BH/include/boost/config.hpp:39:
## /Library/Frameworks/R.framework/Versions/3.3/Resources/library/BH/include/boost/config/compiler/clang.hpp:196:11: warning: 'BOOST_NO_CXX11_RVALUE_REFERENCES' macro redefined [-Wmacro-redefined]
## # define BOOST_NO_CXX11_RVALUE_REFERENCES
## ^
## <command line>:6:9: note: previous definition is here
## #define BOOST_NO_CXX11_RVALUE_REFERENCES 1
## ^
## 1 warning generated.
##
## SAMPLING FOR MODEL 'surv ~ dbinom(density, p)' NOW (CHAIN 1).
##
## Chain 1, Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1, Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1, Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1, Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1, Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1, Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1, Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1, Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1, Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1, Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1, Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1, Iteration: 2000 / 2000 [100%] (Sampling)
## Elapsed Time: 0.177463 seconds (Warm-up)
## 0.105072 seconds (Sampling)
## 0.282535 seconds (Total)## The following numerical problems occured the indicated number of times on chain 1## count
## Exception thrown at line 17: normal_log: Scale parameter is 0, but must be > 0! 1## When a numerical problem occurs, the Hamiltonian proposal gets rejected.## See http://mc-stan.org/misc/warnings.html#exception-hamiltonian-proposal-rejected## If the number in the 'count' column is small, do not ask about this message on stan-users.##
## SAMPLING FOR MODEL 'surv ~ dbinom(density, p)' NOW (CHAIN 2).
##
## Chain 2, Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2, Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2, Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2, Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2, Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2, Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2, Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2, Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2, Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2, Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2, Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2, Iteration: 2000 / 2000 [100%] (Sampling)
## Elapsed Time: 0.18435 seconds (Warm-up)
## 0.112896 seconds (Sampling)
## 0.297246 seconds (Total)
##
##
## SAMPLING FOR MODEL 'surv ~ dbinom(density, p)' NOW (CHAIN 3).
##
## Chain 3, Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3, Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3, Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3, Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3, Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3, Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3, Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3, Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3, Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3, Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3, Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3, Iteration: 2000 / 2000 [100%] (Sampling)
## Elapsed Time: 0.174656 seconds (Warm-up)
## 0.103942 seconds (Sampling)
## 0.278598 seconds (Total)
##
##
## SAMPLING FOR MODEL 'surv ~ dbinom(density, p)' NOW (CHAIN 4).
##
## Chain 4, Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4, Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4, Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4, Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4, Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4, Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4, Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4, Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4, Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4, Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4, Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4, Iteration: 2000 / 2000 [100%] (Sampling)
## Elapsed Time: 0.180502 seconds (Warm-up)
## 0.107807 seconds (Sampling)
## 0.288309 seconds (Total)
##
##
## SAMPLING FOR MODEL 'surv ~ dbinom(density, p)' NOW (CHAIN 1).
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
##
## Chain 1, Iteration: 1 / 1 [100%] (Sampling)
## Elapsed Time: 4e-06 seconds (Warm-up)
## 8.1e-05 seconds (Sampling)
## 8.5e-05 seconds (Total)## Computing WAIC## Constructing posterior predictions## [ 400 / 4000 ]
[ 800 / 4000 ]
[ 1200 / 4000 ]
[ 1600 / 4000 ]
[ 2000 / 4000 ]
[ 2400 / 4000 ]
[ 2800 / 4000 ]
[ 3200 / 4000 ]
[ 3600 / 4000 ]
[ 4000 / 4000 ]## Aggregated binomial counts detected. Splitting to 0/1 outcome for WAIC calculation.precis(m12.2.stan, depth = 2)## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## a_tank[1] 2.12 0.87 0.84 3.60 4000 1
## a_tank[2] 3.04 1.10 1.33 4.72 4000 1
## a_tank[3] 1.00 0.67 -0.07 2.02 4000 1
## a_tank[4] 3.04 1.11 1.21 4.66 4000 1
## a_tank[5] 2.13 0.89 0.59 3.44 4000 1
## a_tank[6] 2.11 0.86 0.74 3.40 4000 1
## a_tank[7] 3.03 1.10 1.34 4.70 4000 1
## a_tank[8] 2.12 0.84 0.83 3.46 4000 1
## a_tank[9] -0.19 0.63 -1.20 0.80 4000 1
## a_tank[10] 2.10 0.85 0.72 3.44 4000 1
## a_tank[11] 1.00 0.67 -0.05 2.06 4000 1
## a_tank[12] 0.57 0.61 -0.35 1.58 4000 1
## a_tank[13] 1.00 0.67 0.00 2.13 4000 1
## a_tank[14] 0.19 0.62 -0.78 1.19 4000 1
## a_tank[15] 2.13 0.87 0.78 3.47 4000 1
## a_tank[16] 2.11 0.84 0.74 3.40 4000 1
## a_tank[17] 2.89 0.78 1.68 4.10 4000 1
## a_tank[18] 2.40 0.68 1.25 3.37 4000 1
## a_tank[19] 2.00 0.57 1.05 2.83 4000 1
## a_tank[20] 3.67 1.00 2.15 5.16 4000 1
## a_tank[21] 2.39 0.67 1.35 3.45 4000 1
## a_tank[22] 2.38 0.66 1.40 3.47 4000 1
## a_tank[23] 2.39 0.66 1.34 3.42 4000 1
## a_tank[24] 1.70 0.54 0.79 2.49 4000 1
## a_tank[25] -1.00 0.45 -1.69 -0.27 4000 1
## a_tank[26] 0.15 0.40 -0.47 0.82 4000 1
## a_tank[27] -1.43 0.49 -2.19 -0.66 4000 1
## a_tank[28] -0.47 0.41 -1.08 0.21 4000 1
## a_tank[29] 0.15 0.41 -0.51 0.80 4000 1
## a_tank[30] 1.44 0.49 0.73 2.28 4000 1
## a_tank[31] -0.64 0.42 -1.29 0.06 4000 1
## a_tank[32] -0.31 0.41 -0.91 0.37 4000 1
## a_tank[33] 3.18 0.78 2.01 4.43 4000 1
## a_tank[34] 2.71 0.65 1.73 3.75 4000 1
## a_tank[35] 2.70 0.66 1.67 3.72 4000 1
## a_tank[36] 2.06 0.51 1.23 2.82 4000 1
## a_tank[37] 2.05 0.49 1.28 2.83 4000 1
## a_tank[38] 3.86 0.96 2.29 5.24 4000 1
## a_tank[39] 2.71 0.67 1.65 3.76 4000 1
## a_tank[40] 2.34 0.56 1.43 3.16 4000 1
## a_tank[41] -1.82 0.48 -2.61 -1.11 4000 1
## a_tank[42] -0.56 0.35 -1.16 -0.03 4000 1
## a_tank[43] -0.46 0.34 -0.98 0.08 4000 1
## a_tank[44] -0.33 0.34 -0.86 0.23 4000 1
## a_tank[45] 0.59 0.34 0.04 1.12 4000 1
## a_tank[46] -0.57 0.35 -1.10 0.00 4000 1
## a_tank[47] 2.06 0.52 1.22 2.84 4000 1
## a_tank[48] 0.00 0.34 -0.55 0.52 4000 1
## a 1.30 0.26 0.88 1.69 4000 1
## sigma 1.62 0.22 1.29 1.96 2729 1#refit the model 12.2 with brms, by specifying different priors for specific coefficients
#(1)without interaction
m12.2.brms1 <- brm(propsurv ~ pred + size,
data = d2,
prior= c(
set_prior("normal(0,1)", class="Intercept"),
set_prior("normal(0,1)", class="b"),
set_prior("cauchy(0,1)", class = "sigma")))## Compiling the C++ model## Start samplingsummary(m12.2.brms1)## Family: gaussian (identity)
## Formula: propsurv ~ pred + size
## Data: d2 (Number of observations: 48)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
## WAIC: Not computed
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## Intercept 0.88 0.04 0.79 0.96 4000 1
## predpred -0.40 0.05 -0.50 -0.30 4000 1
## sizesmall 0.09 0.05 -0.01 0.19 4000 1
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sigma 0.17 0.02 0.14 0.22 4000 1
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).plot(m12.2.brms1,ask=FALSE)
#(2)with interaction
m12.2.brms2 <- brm(propsurv ~ pred * size,
data = d2,
prior= c(
set_prior("normal(0,1)", class="Intercept"),
set_prior("normal(0,1)", class="b"),
set_prior("cauchy(0,1)", class = "sigma")))## Compiling the C++ model
## Start samplingsummary(m12.2.brms2)## Family: gaussian (identity)
## Formula: propsurv ~ pred * size
## Data: d2 (Number of observations: 48)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
## WAIC: Not computed
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## Intercept 0.92 0.05 0.83 1.02 2459 1
## predpred -0.50 0.07 -0.63 -0.37 2435 1
## sizesmall 0.00 0.07 -0.14 0.13 2264 1
## predpred:sizesmall 0.19 0.10 0.00 0.38 2175 1
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sigma 0.17 0.02 0.14 0.21 3294 1
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).plot(m12.2.brms2,ask=FALSE)
#comparing models
loo(m12.1.brms1,m12.1.brms2,m12.2.brms1,m12.2.brms2)## Warning: Some Pareto k diagnostic values are slightly high. See
## help('pareto-k-diagnostic') for details.
## Warning: Some Pareto k diagnostic values are slightly high. See
## help('pareto-k-diagnostic') for details.## LOOIC SE
## m12.1.brms1 -28.91 10.61
## m12.1.brms2 -30.49 12.27
## m12.2.brms1 -28.84 10.68
## m12.2.brms2 -30.48 12.24
## m12.1.brms1 - m12.1.brms2 1.58 4.69
## m12.1.brms1 - m12.2.brms1 -0.07 0.10
## m12.1.brms1 - m12.2.brms2 1.57 4.50
## m12.1.brms2 - m12.2.brms1 -1.65 4.74
## m12.1.brms2 - m12.2.brms2 -0.01 0.23
## m12.2.brms1 - m12.2.brms2 1.65 4.55Q3:For the tomato data set, use both brms and map2stan to fit a model for intleng as a function of species, trt and their interaction, and include shelf as a random effect.
d3 <- read.csv("/Users/xyyue/Documents/RClub/Rclub-rethinking_Xiaoyan.Yue/Assignment_Chapter_09/TomatoR2CSHL.csv")
head(d3)## shelf flat col row acs trt days date hyp int1 int2 int3 int4
## 1 Z 1 B 1 LA2580 H 28 5/5/08 19.46 2.37 1.59 1.87 0.51
## 2 Z 1 C 1 LA1305 H 28 5/5/08 31.28 3.34 0.01 9.19 1.62
## 3 Z 1 D 1 LA1973 H 28 5/5/08 56.65 8.43 2.39 6.70 3.69
## 4 Z 1 E 1 LA2748 H 28 5/5/08 35.18 0.56 0.00 1.60 0.61
## 5 Z 1 F 1 LA2931 H 28 5/5/08 35.32 0.82 0.02 1.49 0.46
## 6 Z 1 G 1 LA1317 H 28 5/5/08 28.74 1.07 6.69 5.72 4.76
## intleng totleng petleng leafleng leafwid leafnum ndvi lat lon
## 1 6.34 25.80 15.78 30.53 34.44 5 111 -9.5167 -78.0083
## 2 14.16 45.44 12.36 22.93 13.99 4 120 -13.3833 -75.3583
## 3 21.21 77.86 13.05 46.71 43.78 5 110 -16.2333 -71.7000
## 4 2.77 37.95 8.08 26.82 33.28 5 105 -20.4833 -69.9833
## 5 2.79 38.11 7.68 22.40 23.61 5 106 -20.9167 -69.0667
## 6 18.24 46.98 23.66 42.35 42.35 5 132 -13.4167 -73.8417
## alt species who
## 1 740 S. pennellii Dan
## 2 3360 S. peruvianum Dan
## 3 2585 S. peruvianum Dan
## 4 1020 S. chilense Dan
## 5 2460 S. chilense Dan
## 6 2000 S. chmielewskii Dan#fit a model with stan for intleng as a function of species, trt and their interaction
#categorical variable for trt, species and shelf
d3$trt2 <- as.numeric(d3$trt)-1 #0 = H, 1=L
d3$species_id <- coerce_index(d3$species)
d3$shelf_id <- coerce_index(d3$shelf)
d3.subset <- d3[,c("shelf_id","trt2","intleng","species_id")]
head(d3.subset)## shelf_id trt2 intleng species_id
## 1 6 0 6.34 4
## 2 6 0 14.16 5
## 3 6 0 21.21 5
## 4 6 0 2.77 1
## 5 6 0 2.79 1
## 6 6 0 18.24 2summary(d3.subset)## shelf_id trt2 intleng species_id
## Min. :1.000 Min. :0.0000 Min. : 0.000 Min. :1.000
## 1st Qu.:2.000 1st Qu.:0.0000 1st Qu.: 9.637 1st Qu.:2.000
## Median :3.000 Median :1.0000 Median :17.255 Median :3.000
## Mean :3.512 Mean :0.5089 Mean :20.340 Mean :2.927
## 3rd Qu.:5.000 3rd Qu.:1.0000 3rd Qu.:28.145 3rd Qu.:4.000
## Max. :6.000 Max. :1.0000 Max. :92.420 Max. :5.000#fit the model
mQ3.stan <- map2stan(
alist(
intleng ~ dnorm(mu,sigma),
mu <- a + a_shelf[shelf_id] + bs[species_id] + bt*trt2 + bst[species_id]*trt2,
a ~ dnorm(20,100),
a_shelf[shelf_id] ~ dnorm(0,sigma_shelf),
bs[species_id] ~ dnorm(0,10),
bt ~ dnorm(0,20),
bst[species_id] ~ dnorm(0,10),
sigma ~ dcauchy(0,1),
sigma_shelf ~ dexp(1)
),
data = d3.subset,
iter = 4000,
warmup = 1000,
chains=4,
cores=2
)## In file included from file10a42593103f.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math/rev/core.hpp:12:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math/rev/core/gevv_vvv_vari.hpp:5:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math/rev/core/var.hpp:7:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/BH/include/boost/math/tools/config.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/BH/include/boost/config.hpp:39:
## /Library/Frameworks/R.framework/Versions/3.3/Resources/library/BH/include/boost/config/compiler/clang.hpp:196:11: warning: 'BOOST_NO_CXX11_RVALUE_REFERENCES' macro redefined [-Wmacro-redefined]
## # define BOOST_NO_CXX11_RVALUE_REFERENCES
## ^
## <command line>:6:9: note: previous definition is here
## #define BOOST_NO_CXX11_RVALUE_REFERENCES 1
## ^
## 1 warning generated.## Warning: There were 1 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup## Warning: Examine the pairs() plot to diagnose sampling problems##
## SAMPLING FOR MODEL 'intleng ~ dnorm(mu, sigma)' NOW (CHAIN 1).
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
##
## Chain 1, Iteration: 1 / 1 [100%] (Sampling)
## Elapsed Time: 4e-06 seconds (Warm-up)
## 0.000257 seconds (Sampling)
## 0.000261 seconds (Total)## Computing WAIC## Constructing posterior predictions## [ 1200 / 12000 ]
[ 2400 / 12000 ]
[ 3600 / 12000 ]
[ 4800 / 12000 ]
[ 6000 / 12000 ]
[ 7200 / 12000 ]
[ 8400 / 12000 ]
[ 9600 / 12000 ]
[ 10800 / 12000 ]
[ 12000 / 12000 ]## Warning in map2stan(alist(intleng ~ dnorm(mu, sigma), mu <- a + a_shelf[shelf_id] + : There were 1 divergent iterations during sampling.
## Check the chains (trace plots, n_eff, Rhat) carefully to ensure they are valid.precis(mQ3.stan)## Warning in precis(mQ3.stan): There were 1 divergent iterations during sampling.
## Check the chains (trace plots, n_eff, Rhat) carefully to ensure they are valid.## 16 vector or matrix parameters omitted in display. Use depth=2 to show them.## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## a 10.67 4.62 2.93 17.67 3228 1
## bt 17.03 4.64 9.27 24.07 3461 1
## sigma 10.61 0.24 10.20 10.97 10429 1
## sigma_shelf 1.42 0.68 0.39 2.44 2920 1pairs(mQ3.stan)
#fit a model with brm
#(1)without interaction
mQ3.brms1 <- brm(intleng ~ 0 + species + trt + (1|shelf),
data = d3,
prior = c(
set_prior("normal(0,10)",class="b"), # sets prior for all b coefficients
set_prior("normal(0,5)",class="b", coef = "trtL"), #set prior for "trtL"
set_prior("cauchy(0,1)", class = "sigma"), #half cauchy prior for sigma
set_prior("normal(0,1)", class = "sd", group = "shelf") #prior for shelf
)
)## Compiling the C++ model## Start samplingsummary(mQ3.brms1)## Family: gaussian (identity)
## Formula: intleng ~ 0 + species + trt + (1 | shelf)
## Data: d3 (Number of observations: 1008)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
## WAIC: Not computed
##
## Group-Level Effects:
## ~shelf (Number of levels: 6)
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sd(Intercept) 1.29 0.49 0.44 2.36 1460 1
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## speciesS.chilense 7.82 1.11 5.57 9.98 1561 1
## speciesS.chmielewskii 11.89 1.10 9.68 14.07 1566 1
## speciesS.habrochaites 15.11 1.11 13.00 17.38 1579 1
## speciesS.pennellii 5.49 1.24 3.13 7.93 1921 1
## speciesS.peruvianum 13.75 1.11 11.54 15.94 1487 1
## trtL 17.13 1.30 14.30 19.50 1302 1
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sigma 10.67 0.24 10.21 11.13 4000 1
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).plot(mQ3.brms1)
#(2)with interaction
mQ3.brms2 <- brm(intleng ~ 0 + trt * species + (1|shelf),
data = d3,
prior = c(
set_prior("normal(0,10)",class="b"), # sets prior for all b coefficients
set_prior("normal(0,5)",class="b", coef = "trtL"), #set prior for "trtL"
set_prior("cauchy(0,1)", class = "sigma"), #half cauchy prior for sigma
set_prior("normal(0,1)", class = "sd", group = "shelf") #prior for shelf
)
)## Compiling the C++ model
## Start sampling## Warning: There were 2 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## Warning: Examine the pairs() plot to diagnose sampling problemssummary(mQ3.brms2)## Family: gaussian (identity)
## Formula: intleng ~ 0 + trt * species + (1 | shelf)
## Data: d3 (Number of observations: 1008)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
## WAIC: Not computed
##
## Group-Level Effects:
## ~shelf (Number of levels: 6)
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sd(Intercept) 1.39 0.56 0.39 2.62 1270 1
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample
## trtH 5.87 1.30 3.24 8.39 1745
## trtL 25.57 1.41 22.61 28.13 1077
## speciesS.chmielewskii 7.75 1.41 4.93 10.45 2120
## speciesS.habrochaites 9.30 1.37 6.59 12.09 1925
## speciesS.pennellii -0.48 1.67 -3.80 2.74 2110
## speciesS.peruvianum 7.41 1.38 4.67 10.01 1817
## trtL:speciesS.chmielewskii -6.26 1.94 -10.05 -2.47 2141
## trtL:speciesS.habrochaites -3.04 1.92 -6.81 0.64 2114
## trtL:speciesS.pennellii -2.63 2.25 -7.09 1.73 1869
## trtL:speciesS.peruvianum -1.97 1.97 -5.85 1.91 2018
## Rhat
## trtH 1
## trtL 1
## speciesS.chmielewskii 1
## speciesS.habrochaites 1
## speciesS.pennellii 1
## speciesS.peruvianum 1
## trtL:speciesS.chmielewskii 1
## trtL:speciesS.habrochaites 1
## trtL:speciesS.pennellii 1
## trtL:speciesS.peruvianum 1
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sigma 10.61 0.24 10.16 11.07 4000 1
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).plot(mQ3.brms2)
#compare the model
loo(mQ3.brms1, mQ3.brms2)## LOOIC SE
## mQ3.brms1 7645.14 64.73
## mQ3.brms2 7638.18 65.98
## mQ3.brms1 - mQ3.brms2 6.96 5.60